% 读取数据文件并转换为数组
data = readtable('加噪声加速度set4');  % 读取表格数据
data = table2array(data);  % 将表格转换为数组

% 采样频率
Fs = 500;  % 采样频率为 500 Hz

% 小波包分解参数
wavelet = 'db4';  % 小波基函数
level = 3;        % 分解层数
cutoff_frequency = 120;  % 截止频率 120 Hz
data = lowpass(data, cutoff_frequency, 2000);

% -----------------------------
% 不加白噪声的小波包分解
% -----------------------------
wpt_no_noise = wpdec(data, level, wavelet);  % 不加噪声的小波包分解

% 获取第四层的节点
tn_no_noise = leaves(wpt_no_noise, level); % 获取第四层节点

% 初始化存储不加噪声情况下的能量熵
energy_entropy_no_noise = zeros(1, length(tn_no_noise));

% 计算总信号能量
total_energy = sum(data.^2);

% 计算不加噪声的每个节点的能量熵
for i = 1:length(tn_no_noise)
    node_signal = wprcoef(wpt_no_noise, tn_no_noise(i));
    
    % 计算节点能量
    node_energy = sum(node_signal.^2);
    
    % 计算能量熵
    energy_ratio = node_energy / total_energy;
    energy_entropy_no_noise(i) = -energy_ratio * log2(energy_ratio + eps);
end

% -----------------------------
% 加白噪声的小波包分解
% -----------------------------
% 添加白噪声
noise_level = 0.3;
data_noisy = data + noise_level * randn(size(data)); % 加噪声的信号
total_energy2 = sum(data_noisy.^2);
% 小波包分解加噪声信号
wpt_with_noise = wpdec(data_noisy, level, wavelet);
tn_with_noise = leaves(wpt_with_noise, level); % 获取第三层节点

% 初始化存储加噪声情况下的能量熵
energy_entropy_with_noise = zeros(1, length(tn_with_noise));

% 计算加噪声的每个节点的能量熵
for i = 1:length(tn_with_noise)
    node_signal = wprcoef(wpt_with_noise, tn_with_noise(i));
    
    % 计算节点能量
    node_energy = sum(node_signal.^2);
    
    % 计算能量熵
    energy_ratio = node_energy / total_energy2;
    energy_entropy_with_noise(i) = -energy_ratio * log2(energy_ratio + eps);
end

% -----------------------------
% 使用 K-means 聚类分类信号类型 (不加噪声情况)
% -----------------------------
num_clusters = 3;  % 我们将能量熵分为3类：低频信号，正常信号，噪声信号

% 执行 K-means 聚类
[idx_no_noise, C_no_noise] = kmeans(energy_entropy_no_noise', num_clusters);

% 初始化分类数组
signal_types_no_noise = strings(1, length(energy_entropy_no_noise));

% 将每个聚类分组映射为信号类型
for i = 1:length(idx_no_noise)
    if idx_no_noise(i) == find(C_no_noise == min(C_no_noise))  % 最小类中心对应低频信号
        signal_types_no_noise(i) = '低频信号';
    elseif idx_no_noise(i) == find(C_no_noise == max(C_no_noise))  % 最大类中心对应噪声信号
        signal_types_no_noise(i) = '噪声信号';
    else  % 中间类中心对应正常信号
        signal_types_no_noise(i) = '正常信号';
    end
end

% -----------------------------
% 使用 K-means 聚类分类信号类型 (加噪声情况)
% -----------------------------
% 执行 K-means 聚类
[idx_with_noise, C_with_noise] = kmeans(energy_entropy_with_noise', num_clusters);

% 初始化分类数组
signal_types_with_noise = strings(1, length(energy_entropy_with_noise));

% 将每个聚类分组映射为信号类型
for i = 1:length(idx_with_noise)
    if idx_with_noise(i) == find(C_with_noise == min(C_with_noise))  % 最小类中心对应低频信号
        signal_types_with_noise(i) = '低频信号';
    elseif idx_with_noise(i) == find(C_with_noise == max(C_with_noise))  % 最大类中心对应噪声信号
        signal_types_with_noise(i) = '噪声信号';
    else  % 中间类中心对应正常信号
        signal_types_with_noise(i) = '正常信号';
    end
end

% -----------------------------
% 绘制不加噪声与加噪声情况下的能量熵对比
% -----------------------------
figure;

% 不加噪声的能量熵
subplot(2,1,1);
plot(0:length(energy_entropy_no_noise)-1, energy_entropy_no_noise, '-o', 'LineWidth', 2, 'MarkerSize', 8);
hold on;

% 分类标记
for i = 1:length(tn_no_noise)
    text(i-1, energy_entropy_no_noise(i) + 0.0002, signal_types_no_noise(i), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
end
title('不加噪声的第三层节点能量熵及信号类型 (K-means 聚类)');
xlabel('节点编号');
ylabel('能量熵');
grid on;

% 加噪声的能量熵
subplot(2,1,2);
plot(0:length(energy_entropy_with_noise)-1, energy_entropy_with_noise, '-o', 'LineWidth', 2, 'MarkerSize', 8);
hold on;

% 分类标记
for i = 1:length(tn_with_noise)
    text(i-1, energy_entropy_with_noise(i) + 0.0002, signal_types_with_noise(i), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
end
title('加噪声的第三层节点能量熵及信号类型 (K-means 聚类)');
xlabel('节点编号');
ylabel('能量熵');
grid on;
